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ABSTRACT 


We use a block Lanczos algorithm for computing a few of the smallest 
eigenvalues and the corresponding eigenvectors of a large symmetric matrix rather 
than computing all the eigenvalue-eigenvector pairs. The basic Lanczos algorithm 
generates a similar matrix which is block tridiagonal from a given large symmetric 
matrix. The size of the generated tridiagonal matrix depends upon the number of 
the smallest eigenvalues to be computed. The result is savings in computations and 
storage. The block Lanczos algorithm is well-suited for problems involving multiple 
eigenvalues. 

In this thesis, we develop the block Lanczos algorithm to estimate the 
direction-of-arrival (DOA) of a point source based on the observations measured at 
a linear array of sensors and compare the performance with that of a single vector 
Lanczos algorithm. The results of the computer simulation experiments conducted 


with this method are presented and discussed. 
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I. INTRODUCTION 


It is necessary to estimate the power spectral density (PSD) function in 
several signal processing problems. Of particular interest is the estimation of 
spectral lines in noise. One related problem is the direction—of—arrival estimation of 
point sources based on array measurements. Methods based on subspace estimation 
have been proposed in this regard [Ref. 9]. Much of the basic literature on these 
methods has been borrowed from functional and numerical analysis [Refs. 8, 21]. 
Determination of the subspaces needs a complete eigendecomposition of the given 
data autocorrelation matrix. The general eigendecomposition of an nxn matrix 
requires O(n?) computations. Often we only need to compute a few of either the 
smallest or the largest eigenvalues and the corresponding eigenvectors of a large 
symmetric autocorrelation matrix rather than all the eigenpairs of the matrix. In 
this research, we develop an algorithm for computing a few of the smallest 
eigenvalues and the corresponding eigenvectors and apply this to high resolution 
spectral estimations problems. 

The algorithm we develop is an extension of the method of minimized 
iterations due to Lanczos [Ref. 9]. Paige {Ref. 4] experimented with the Lanczos 
algorithm and found that a few of the extreme eigenvalues of a tridiagonal matrix 
would often converge rapidly to the similar eigenvalues of a real symmetric matrix 
R much before the entire set of eigenvalues are computed. 

The block Lanczos algorithm is an extension of the basic single vector Lanczos 
algorithm. We iterate with a block of vectors rather than with a single vector, and 


generate a reduced block tridiagonal matrix that is similar to the original 


autocorrelation matrix R. The basic idea is that the eigenvalues of the block 
tridiagonal matrix are approximately the same as the extreme eigenvalues of R. 
Several researchers have reported the block Lanczos algorithm and its variants, in 
particular, Cullum and Willoughby [Ref. 1], Kahan and Parlett [Ref. 6] and Golub 
and Underwood [Ref. 4]. 

The objective of this thesis is to develop a block Lanczos algorithm to extract 
a few of the smallest eigenvalues and then estimate the corresponding eigenvectors. 
The smallest eigenvalues are said to correspond to the noise subspace of the spectral 
Or array measurements. The proposed algorithm will be used to estimate the 
spectral lines which may represent the direction-of-arrival of point sources in low 
signal-to-noise ratio environments. The block Lanczos algorithm will have to be 
compared with the single vector case with respect to the spectral line estimation 
performance. With these goals in mind, we now proceed to discuss how the thesis 
is organized. 

In Chapter II we introduce and summarize the single vector Lanczos 
algorithm. Complete and selective reorthogonalization of Lanczos vectors is 
discussed there. In Chapter III, we describe and develop the block Lanczos 
algorithm and present the results of the experiments carried out with a computer 
program implementing this algorithm. Also, in this chapter, we compare the 
performance of DOA estimation of the block Lanczos algorithm with that of the 
single vector Lanczos algorithm. Finally, in the last chapter, we discuss and 
summarize the results of the Lanczos method and also make some recommendations 


for future work. 


H. SINGLE VECTOR LANCZOS ALGORITHM 


The single vector Lanczos algorithm is used to tridiagonalize a real symmetric 
matrix R. This algorithm is based on the concepts such as Krylov sequences and 
subspaces, orthogonal projections of matrices, and Ritz vectors [Ref. 1]. In this 
chapter we describe the general Lanczos recursion and a practical Lanczos 
algorithm. Issues related to the reorthogonalization of Lanczos vectors are also 


discussed. 


A. LANCZOS RECURSION 
Given an nxn real symmetrix matrix R and an arbitrary nx] unit norm vector 
k,, we can obtain a sequence of vectors defining the n dimensional Krylov subspace 


as follows [Ref. 1] 


k, = hie 
kote NK 
مر دبا‎ ks (2-1) 


We can now form the Krylov matrix of rank m as follows 


K,, = [ky k, kg + + + Kp] 
RR Tk (2-2) 


and the Krylov subspaces J£", m=1,2,---,n, are then defined as 
A™(R,k,) = span {ki Rk Rk Rk (2—3) 


We now proceed to obtain a tridiagonal matrix which has the same extreme 
eigenvalues as the given symmetric matrix R. If the columns of a matrix Q,, are 


orthonormal, then the matrix 
T 
E = Q, RQ, (2-4) 


is an mxm tridiagonal matrix which is an orthogonal projection of R onto the space 
spanned by the columns of Q,. Lanczos [Refs. 1, 9, 14] proposed a recursion to 
generate the columns of the matrix Q, — [q, 4, *** q,]. The vectors q,, q2, °°°, Un 
are referred to as Lanczos vectors. Thus, for a given real symmetric matrix R, the 
Lanczos recursion produces a tridiagonal matrix T,, as a projection of R onto the 
corresponding subspace span{Q,,}, spanned by the Lanczos vectors generated. 
These subspaces are in fact Krylov subspaces. [Ref. 1] 

The classical Lanczos recursion for R starts with an nx1 initial Lanczos vector 


qı which is randomly generated and normalized. Initially we define 0,2 0 and qo= 0. 


We then compute a,, §;,, and the Lanczos vectors q, for i=1,2,- - -,m as follows: 
T 
0; — q,Rq; (2-5) 
يې ز01‎ = Ra; — od; — idi- (2—6) 


T 
bini 7 (Zan) 


The Lanczos matrix T,, is generated by appropriately arranging a, and £;,,, where 
0; are the diagonal elements and 5;,, are the subdiagonal elements of the tridiagonal 


matrix, so that 


TRGI) a for 1<i<m, (2-8) 
and 
1 111711 ال‎ 3711 20707 11 (2-9) 
resulting in 
a, By 
By O û; 
Pz 03 
Ta = (2—10) 
Pm 
Bn n | 


The vectors a; 


q; and /;q;., in Eqn(2—6) are the orthogonal projections of vectors 


Rq; onto the two most recently generated Lanczos vectors q; and q;- Thus, 


1 
updated Lanczos vector q;,, is computed by orthogonalizing the vector Rq; with 


respect to previously computed Lanczos vectors q; and q;.,. The classical Lanczos 


recursion can be condensed into matrix notation by a group of Lanczos vectors Q,, 


and the Lanczos matrix T,. We obtain the following matrix equation 
T 
RQ, = Ql + زې‎ (2-11) 


where e, is the unit vector whose ¿ element is 1 and whose other elements are 0. 
Note that, in this recursion, R is never modified. Also, storage is needed only for 
the Lanczos vectors q;-,, q; and q;,,, the Lanczos matrix T,, and the given matrix 
R. 

The classical Lanczos procedure attempts to maintain the orthogonality of the 
Lanczos vectors. Due to the roundoff and other numerical errors, however, 
orthogonality between the Lanczos vectors can only be maintained by incorporating 
some kind of explicit reorthogonalization as the Lanczos vectors are computed. 
Note that the reorthogonalization of these vectors requires extra storage for keeping 
all of the Lanczos vectors as well as additional computations. We discuss the 


reorthogonalization in a later section. 


B. PRACTICAL LANCZOS ALGORITHM 

In this section we focus on a Lanczos algorithm with no explicit 
reorthogonalization incorporated. The loss of orthogonality, if it occurs, may result 
in Spurious eigenvalues. Although the orthogonality of the Lanczos vectors is lost, 
some of the eigenvalues of R will suill appear as eigenvalues of the Lanczos matrix if 
we make the tridiagonal matrix large enough. Since we are generally interested in 
only a few of the smallest eigenvalues and their corresponding eigenvectors, the loss 
of MINIS does not critically affect finding approximate eigenvalues and 


corresponding eigenvectors of R [Ref. 1]. Nevertheless, it is possible to find accurate 


eigenvalues and eigenvectors by using the method in an iterative scheme using no 
reorthogonalization at all, even in the face of total loss of orthogonality [Ref. 12]. 
Now, we present the practical single vector Lanczos algorithm. 

To generate the tridiagonal Lanczos matrix T,, we should compute a; and (;,, 
which are the diagonal and sub-diagonal elements of the tridiagonal matrix, where 
i=1,2,---,m. Given an nxn real symmetric matrix R, we start with an nx] initial 
arbitrary vector q, which is normalized such that |q,|.=1 as in Eqn(2—1). We now 


define an intermediate vector 


u, = Rq; (2-12) 


and initialize a5—0 and 55-0. The single vector recursion is then carried out for 


i=1,2,---,m and can be summarized as follows: 
T 
a = a (2-13) 


T = +y ww; (2—15) 


qii = W;/T ys (2=16) 
T 

6, = a;Rqj,, (EN 

Uj, = Kqi4 — fq (2-18) 


where u, and w, are the orthonormal projection vectors for the Lanczos vector q;. 
Here Eqn(2—13) and Eqn(2—17) use the modified Gram-Schmidt orthogonalization 
to compute the coefficients o; and fi,,, respectively. If q; is orthogonal to q;_,, 
then q;,, will theoretically be orthogonal to both q;., and q;. This algorithm is 
quite popular despite the requirement for small amounts of extra computations 
[Ref. 12]. Now the tridiagonal matrix T,, is generated by simply filling it with a, 


and Û, for its entries as shown in Eqn(2—10). 


1. Eigenvalue Computation 
To find the eigenvalues y, of the mxm Lanczos matrix T,,, we may use 
the bisection method and Sturm sequencing [Ref. 5]. Actually, one could obtain 
both eigenvalues and eigenvectors of T,, by using such methods as the QR algorithm 
[Ref. 16]. However, since we need only a few of the smallest eigenvalues of Ta, we 
choose the bisection method to compute them as detailed in the following. 


Given the tridiagonal matrix T,, we define the characteristic 


polynomials po(u), p,(u), * * ^, pg (u) as 
po(u) = 1 (2-19) 
pj(4) ^ det(T, — Hl), (2-20) 


for j=1,2,---,m. Expanding the determinant in Eqn(2—20) yields the recursive 


expression [Ref. 5: pp. 305—307]: 


p;(u) 7 (o5 — )pj- (1) — Bipj-o(4), — j92,3,- m. (2-21) 


The zeros of the polynomial p,(j) are the eigenvalues of the tridiagonal matrix Tm 
Here we are interested in only a few of the smallest eigenvalues of R. In order to 
compute these values we first define a range [x, y] in which all the desired 
eigenvalues lie. We then carry’ out the following iteration to implement the Sturm 


sequencing property. If p(x)p(y) < 0 and x < y, then the iteration 


IO 


Inm dy 
[iss 5 
y=p if p,(x)p, (y) « O (222) 


x=p if pg(x)pg(y) 2 0 


is guaranteed to converge to a zero of p, (1), i.e., to an eigenvalue of Ta. The value 


e is the machine unit roundoff error and the limits on the range [x, y] are given by 


o nn A 
(2-23) 
y = max a; + [8] + |8;.,| 
1 


where we have by = §,,,, = 0. The bisection method computes the eigenvalues with 
small relative error, regardless of their magnitude. This is in contrast to the 
tridiagonal QR iteration, where the computed eigenvalues can have only small 


absolute error. 


2. | Eigenvector Computation 
Now, we need to compute the eigenvectors of R to complete the 
eigenpair estimation problem. There are two techniques to compute the 
approximate eigenvectors of R. The first technique uses the Ritz vector X; and the 
corresponding algorithm uses the relationship [Ref. 14] 


ا 


CA Po: (2-24) 


where Z; is one of the eigenvectors of T,,, Q,, is the matrix of Lanczos vectors, and 


s < 


4 - 


The second technique to obtain the eigenvectors of R corresponding to 
the selected eigenvalues of T, uses the Rayleigh quotient iteration. The iteration 
makes the Rayleigh quotient r(x,) converge to u; which is one of the smallest 
eigenvalues of T,. The Rayleigh quotient of an eigenvector x, which minimizes 


|(R — r,I)x, |; is given by 


1 
x, Rx, 


1 


(2-25) 
X k Xk 


ry, — r(x,) — 


As r(x,) approaches an eigenvalue u; of Tw then the solution to (R — r,I)x, — b, 
will be an approximate eigenvector by using the inverse iteration theory where b, is 
a vector close to zero [Ref. 14]. 

We now briefly present the inverse iteration algorithm. First, we pick 


an arbitrary unit vector xy; then the iteration proceeds as follows 
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rD = xD TRO, ey (i) 1 
(2-26) 
solve (Ro OL Y 
Xj = Ye Yel 
If r{) converges to one of the smallest eigenvalues y; of T,, where i=1,2,---,s, then 


we stop the iteration and the corresponding vector x{1) is the required eigenvector 
of R. 

Having discussed the Lanczos algorithm and the methods to compute the 
desired eigenpairs, we now proceed to address some issues related to the 


reorthogonalization of the Lanczos vectors. 


C. REORTHOGONALIZATION 
As mentioned in the previous sections, the Lanczos vectors q,, Qo, °°°, Qm lose 
mutual orthogonality as the number of steps m in the algorithm increases. The 


requirement that ala, = I, is then destroyed by the roundoff errors and the 


m 
algorithm is described as unstable. A few steps later, the matrix of Lanczos vectors 
Q,, may not even have full rank, t.e., the Lanczos vectors may become linearly 
dependent. As a result there is no guarantee that T,, will bear any useful 
relationship to R. And the orthogonality among vectors q,, qo, * : *, q, disappears. 


The ideal Lanczos algorithm should terminate (f, = 0) for some m € n, but in 


practice the process goes on forever computing more and more spurious eigenvectors 


Io 


for each correct eigenpair it discovers. [Ref. 14] 

Thus, sometimes we need to maintain the orthogonality of the Lanczos vectors 
for finding more accurate eigenvalues and eigenvectors of R and to avoid the 
computation of any spurious pairs. We now add the reorthogonalization step in the 
single vector Lanczos procedure presented in Section A and B of this Chater. There 
are two techiques for reorthogonalization, namely, complete reorthogonalization 


and selective reorthogonalization. These are discussed in the following. 


1. Complete Reorthogonalization 
Complete  reorthogonalization incorporates a Householder matrix 
computation into the Lanczos algorithm for producing Lanczos vectors that are 
orthogonal to within the working accuracy. This is effective at maintaining the 
stability of the system. The following step is inserted into the Lanczos algorithm 


(Eqns(2—13)-(2—18)) after computing a projection vector w;: 
W; — W; — q;(q; wj), ز10‎ 1 2| (2-27) 


thus w; is explicity orthogonalized against q; and q;.,. If a Householder matrix P; is 
determined so that [Po P -++ Pi] [wp wp +++, w;] is upper triangular, then it 
follows that the (¿+1)st column of [P, +++ P;] is the desired unit vector. An 
example of a complete reorthogonalization Lanczos scheme is summarized below. 
First, we determine the initial Householder matrix Pg = 1 — O so that 
P¿w, = e From Eqn(2—5) we can then obtain a, = qRa, and implement the 


following recursion for i=1,2,---,m—1: [Ref. 5: pp. 334-345] 
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r; = (R— 0,1)q; — B;-14;-1  (B040=0) 


(2-28) 


= 
| 
5 
"d 
E 
= 


: T; T 
determine P; = I — 2v;vį/v;V; 


T 
Pow; = (Wı کس از‎ Û, 0, E 0) 


Gee ge Be 
ji 
Qi. — q;+1RQ;+1- 


In Eqn(2—28) the computed Lanczos vectors q; are now mutually orthogonal to the 
working precision which follows from the roundoff properties of Householder 
matrices [Ref. 5]. However, the complete reorthogonalization of Lanczos vectors 
requires extra computations and all the computed vectors q) qə *** Q;, W, Wo *** W; 


1 


need to be stored. It negates any advantage of the Lanczos algorithm. 


2. selective Reorthogonalization 
Selective reorthogonalization is computationally more efficient. A small 
modification to the exact version of the Lanczos algorithm ensures that the Lanczos 
vectors qp Q5, q4:-- maintain the orthogonality. At the i step the Lanczos 
algorithm (Eqns(2—13) — (2—18)) produces the matrix of Lanczos vectors Q;, the 
tridiagonal matrix T;— QRQ,, and the residual vector w; — (RQ; — Q;T;)e;. The 


1 


equation to compute the residual vectors RQ, = Q;T, + wie; can then be rewritten 
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as 
T 


where F; accounts for the roundoff errors. When the Lanczos vectors lose their 


orthogonality we have 
T 
l1- QiQ;l < k; (2-30) 


where it is required to keep k, < k for some k in the interval (ne, 0.01) and e is the 
relative precision of the arithmetic. 
The loss of orthogonality goes hand in hand with the convergence of a 
Ritz pair. Suppose that the symmetric QR algorithm [Ref. 16] is applied to T; and 
renders the computed Ritz values @,, ---, 0; and a nearly orthogonal matrix of 
eigenvectors Z,=[Z,,Z),--:,Z,| of T,. Then the Ritz vectors which are the 
approximate eigenvectors of R are given by 
X= xo EE O (2—31) 


171 


We can show that the absolute value of the inner product between q;,, and x; is 


c [R]; 
T 
Iaiaxl *—pETTZZAT- (2—32) 
1 ij 


where the denominator term can be approximated as 
۱/١ |z;] S |Rx; — 0;x;1, (2-33) 
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for j=1, 2, ---, i. The recently updated Lanczos vector q;,, is forced to have a 
unwanted component in the direction of any converged Ritz vector. Consequently, 
instead of orthogonalizing q,,, against all of the previous Lanczos vectors, we can 
achieve the same effect by orthogonalizing it against the converged Ritz vectors 
[Ref. 5]. A selective reorthogonalization method based on this technique is 
discussed in Parlett and Scott [Ref. 15]. In this method, a computed Ritz pair (@,x) 


is considered a good approximation if it satisfies 


[Rx — 0x1, 2 V e| RI (2-34) 


where ¢ is the machine precision constant. After computing q;,,, it is 
orthogonalized against each good Ritz vector [Ref. 5]. Selective reorthogonalization 
prevents the computation of many spurious eigenvectors. This means that the extra 
computations and storage required in this method are less than those required to 
implement the complete reorthogonalization since there are usually fewer Ritz 
vectors than Lanczos vectors. When all of the eigenvalues of R are required, 
however, then the selective reorthogonalization procedures are too expensive to 


implement [Ref. 12]. 


D. RESULTS 
In this section we will discuss the experimental results based on computer 
simulation of the single vector Lanczos algorithm for estimating the 


direction-of-arrival of point sources. 
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1. | Experimental Set Up 
Consider that we receive signals at a linear array containing M equally- 
spaced sensors. ‘The signal is modeled as a sum of / sinusoids, the normalizing 
spatial frequencies of which are propotional to their bearings from 0.0 (= 0°) to 
0.5 (= 90°), in additive random noise with fixed variance. The received signal at 


the mt! sensor location is then given by 


l 


2(m) =} A;cos(27 fm) + n(m) (2-35) 


where A; is the amplitude of the iP sinusoid, f, is the normalized spatial frequency 
of the i" signal which represents the bearing, and n(m) is zero-mean white 


2 


Gaussian noise with variance c^. The amplitude of the signal A; and the noise 


variance o? will determine the signal to noise ratio of the i" signal, 


Ai 





SNR; = 10 log( — j (2-36) 


Figure 1 shows the block diagram of the direction-of-arrival estimation 
algorithm considered in this work. After measuring the signals received at each 
sensor, the sensor output is passed through a bank of bandpass filters. This helps 
prefilter the noise over the selected frequency band of frequencies and provides some 
initial processing gain [Ref. 22]. The autocorrelation matrix is then computed by 
taking data from the outputs of corresponding bandpass filters at each sensor. An 


nxn autocorrelation matrix is generated as follows 
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M-1-k 
Ra(k) == ) z(m-k)z(m) fortück&n-l. (2-37) 


m=0 


We now have computed an nxn autocorrelation matrix of z(m), R,,, by 
using M data samples. The eigenvectors x, of R,, corresponding to the smallest 
eigenvalues are computed by using the single vector Lanczos algorithm. The power 


spectral density estimates are computed as: 


SU TEN Ncc acct (2—38) 


where x;; are the elements of the i^ eigenvector x;. 
2. DOA Estimation 

A few of the extreme eigenvalues and their corresponding eigenvectors of 
a large real symmetric matrix R can be obtained by using the single vector Lanczos 
algorithm. Those eigenpairs are the approximate eigenpairs of R. Each eigenvector 
has spectral information to determine the bearing of a source relative to an array of 
sensors. The physical implementation of a DOA estimation scheme is shown in 
Figure 1. The correlator generates an autocorrelation matrix of the received signal. 
The Lanczos algorithm and the eigendecomposition produce the noise subspace 
eigenvectors that estimate the spatial power spectral density (PSD) function, given 
in Eqn(2—38). The PSD function represents the direction-of-arrival of point 
sources as spectral peaks. In this thesis we have used three different ways of 


estimating the PSD function: the individual eigenvector spectra, the algebraic 
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averaging of a few eigenvectors, and the spectral multiplication of the individual 
spectra. 

The computer simulation experiment consists of an equally-spaced array 
of 100 sensors arranged in a linear fashion. By using a filter bank as shown in 
Figure 1, the signals have known temporal frequency with unknown bearings. The 
autocorrelation matrix size is chosen to be 25x25. The number of eigenpairs 
computed is chosen according to the number of iterations in the Lanczos algorithm. 
We have used 5 dB, 0 dB, —5 dB and —10 dB SNR cases for the direction-of-arrival 
estimation. The results indicate the ability of the algorithm to determine the 
number of targets and bearing resolution for various directions and different SNRs. 
In each case we have used the five smallest eigenvalues and their corresponding 
eigenvectors for computing the PSD function. (We can, however, choose more 
eigenvectors at the expense of more computations.) 

We can use eigenvector averaging, which is the algebraic averaging of 
the computed eigenvectors corresponding to the smallest eigenvalues of R. It has 
improved the performance compared to that of the individual eigenvector spectra. 
Further improvement in results, however, can be obtained by spectral multiplication 


of the individual eigenvector spectra. The estimate is given by 


— 
STOP E (2-39) 


1=1 


where J is a predetermined number (J < m < n). 
We now consider several examples to study the estimation performance 
of the Lanczos algorithm and the consequent eigenpair computation. Most of the 


results use the Ritz vector method with no reorthogonalization. Example 1 is the 
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detection of a single target at 9° for different SNRs. Figure 2 shows the overlayed 
individual spectra of the five eigenvectors corresponding to the smallest eigenvalues 
at an SNR of 5 dB. Note that the spectrum of each eigenvector has several spurious 
peaks, but each eigenvector has a common peak at the true bearing 9°. Figure 3(a) 
shows a plot of the PSD function which is the average of 5 eigenvectors. The 
averaging improves the estimation performance considerably. However, further 
improvement was obtained by using the spectral multiplication as defined in 
Eqn(2—39) and the result is shown in Figure 3(b). As can be seen, even though it 
has two small spurious peaks, it has greatly improved the nulls and the peak at 9°. 
In the remainder of the thesis, we use the spectral multiplication method to 
compute the PSD function. Figure 4 and 5 show the results at 0 dB and —5 dB SNR 
respectively. More spurious peaks are observed as the SNR is lowered. At an SNR 
of —10 dB (see Figure 6(a)), several spurious peaks are seen which are almost as 
large as the true bearing. Improved performance is obtained as shown in Figure 
6(b) by using more eigenvectors (7 eigenvectors) in this case in the spectral 
multiplication. 

In Example 2 we have 3 targets at 34°, 36° and 54°. Notice that two 
targets are very closly spaced in bearing. Figure 7(a) is obtained by using the Ritz 
vector method and Figure 7(b) is obtained by using the Rayleigh quotient iteration 
to compute the eigenvectors at 0 dB. Although the two targets are very closely 
spaced, good resolution is clearly achieved and almost no spurious peaks are seen in 
both results. Figure 8(a) and Figure 8(b) are obtained by using no 
reorthogonalization and with complete reorthogonalization respectively at an SNR 
of —5 dB. As mentioned earlier, loss of orthogonality does not affect the results for 


a few eigenvectors in the single vector Lanczos procedure. The result of 


s 


reorthogonalization is almost the same as the result of no reorthogonalization. 
Resolution is achieved in both cases at —5 dB, but several spurious peaks are 
present. Figure 9(a) shows the result using 3 eigenvectors while Figure 9(b) 
indicates the performance using 5 eigenvectors at an SNR of —10 dB. Note that 
sufficient spectral resolution is not acheived to discriminate the targets located at 
34° and 36°. A number of spurious peaks are higher than the peak at 36° making it 
impossible to accurately determine the number of targets as well as their locations. 
The results in this chapter indicate that the eigenvectors found using the 
single vector Lanczos algorithm are sufficiently accurate to determine the spectrum. 
The spectral multiplication scheme achieves the best spectral estimation 
performance. This algorithm provides savings in computations and storage because 
it needs to compute only a few of the extreme eigenvalues and eigenvectors of a 
large real symmetric matrix. Loss of orthogonality is not critically affected when we 
need to find only a few eigenvalues and eigenvectors of a large autocorrelation 


Matrix. 
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III. BLOCK LANCZOS ALGORITHM 


In the previous chapter we presented the single vector Lanczos method. The 
single vector method can be used to find a few extreme eigenvalues of any given real 
symmetric matrix R. However, this method does not determine the multiplicities of 
the eigenvalues directly. Besides, it does not determine a complete basis for the 
invariant subspace corresponding to any such multiple eigenvalue. Here we consider 
an alternative approach which directly determines the multiplicities of the 
eigenvalues and the corresponding eigenvectors |Ref. 1]. 

The method known as the block Lanczos algorithm is an extension of the 
Lanczos algorithm in which a block of vectors rather than a single vector is iterated 
[Ref. 4]. We produce a block tridiagonal matrix in place of the usual tridiagonal 
matrix produced in the single vector Lanczos method. The block Lanczos method 
can be used in a manner proposed by Paige [Ref. 10]. That is, one can compute a 
sequence of estimates of the eigenvalues of the matrix R from the block tridiagonal 
matrix. 

several researchers have worked on the block Lanczos algorithm, in particurlar 
Kahan and Parlett [Ref. 6], Cullum and Willoughby [Ref 1], and Golub and 
Underwood [Ref. 4]. In this chapter we describe the basic idea of the block Lanczos 
method, develop the algorithm, present several simulation results, and compare 


them with those of the single vector algorithm. 


Dd 


A. BLOCK LANCZOS METHOD 

Given a large real correlation matrix R of size nxn, we can generate a banded 
tridiagonal matrix T, of size qsxqs, where qs<n, q is the block size, and s is the 
number of the Lanczos blocks. Starting from an initial nxq orthonormal matrix Q,, 
the purpose is to compute a sequence of mutually orthonormal nxq matrices Q,, Qs, 
Q4, -+ +, Q, such that the space of vectors spanned by the columns of these matrices 
contains the columns of the matrices Q,, RQ,, R?Q,, ..., R*Q!Q,, where 0«q«2 and 


1<s<q. Note that usually we have ps« «n. The block Lanczos algorithm can be 


summarized as follows [Ref. 1]: For i=1,2,---,s, compute 
Ay = Q;(RQ; — Q;1B;) (3-1) 
P; = RQ; -— QA; — Q;-,B; (3-2) 
Q;+¡B¡+, = P; (QR factorization of P;). (3-3) 


In this procedure, Q;,, is orthonormal to all previous Q;. The purpose of the block 
Lanczos algorithm is to find A,, B;,, and Q.,,, where A, and B,,, are the element 
matrices of the desired tridiagonal matrix, and Q;,, is the (i+1)S* orthonormal 
Lanczos block. The blocks Q,, for i=1,2,---,s, form an orthonormal basis of the 


Krylov subspace, defined as 


%*(Q,,R) = span{Q,,RQ,,R7Q,, °°°, R°1Q,} (3-4) 


corresponding to the first block Q,. The matrix Q, defined as a catenation of the 
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Lanczos blocks 


Q = [Q; Q Q3 +++ Q, (35) 


is orthogonal [Ref. 4], :.e., Q'Q = 1 on Q! — Q'. We then generate the banded 


tridiagonal matrix T, using A, and B;: 


A, Bj 
B, A, Bl 0 
B3 Az 
T,- | (3-6) 
0 Bl 
B, ius 


The elements A; and B; are qxq coefficient matrices, so that T, is a qsxqs banded 
matrix. The off-diagonal blocks, B;, are upper triangular matrices and the main 
diagonal blocks, A;, are symmetric matrices so that T, itself is a symmetric matrix. 


Also, the banded tridiagonal matrix T, can be determined by 
T - 
T, = Q RQ. (3—7) 


There are basically two different types of block Lanczos procedures, namely, the 
iterative procedure and the noniterative procedure. The noniterative procedure 
follows along the lines of the single vector Lanczos procedure. In this method a 
sequence of blocks {Q,,Q),۰۰۰} are generated where the length of the sequence can 


be determined by the size of T, and the amount of storage available. The generated 
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Lanczos blocks may or may not be reorthogonalized. The iterative procedure uses 
the block recursion to generate the block tridiagonal matrix. First, the relevant 
eigenvalues and eigenvectors of these block tridiagonal matrices are computed. 
Then the corresponding Ritz vectors are computed and used as updated 
approximations to the desired eigenvectors. If convergence has not been achieved in 
k iterations (see Step 2 below), more iterations are performed these updated 
eigenvectors until the procedure has converged. 


The following steps show the basic iterative block Lanczos procedure: 


Step 1. For k=1 start with an initial arbitrary nxq block Qk where the 
columns of Q* are orthonormal. 

Step 2. — Compute PK = RQ‘ — QKAK using Qk where Ak — (Q5) RQ* and 
use the norms of the columns of P* to check for convergence. If 


convergence has occured, then stop; otherwise, go to Step 3. 


Step 3. Generate a sequence of blocks Qi using the recursion in Eqn(3—2) 
and Eqn(3—3) for j=2, 3, ---, s. Use the coefficient matrices A‘ 
and B¥,, to define the real symmetric block tridiagonal matrix Tf. 

Step 4. Compute the q algebraically smallest eigenvalues of TX and the 
correspoding eigenvectors Y* where Y* = {y}, y§, ‘°, yk}. 


k+ 
Step 5. Obtain the new Lanczos block Q, given by 


k+1 kk 
Q; -QY (38) 
k k k k 
where Q = {Q,, Qo,::-, Q,}. Increment k to k+1 and go to Step 2. 


From the above procedure we can generate the block tridiagonal matrix T, and also 
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compute the eigenvectors Y* which approximate the eigenvectors of R. As noted, 
Step 2 provides a check for the convergence of the procedure. 

If it were the case that qs = n, then T, would be similar to R and the 
eigenvalues of T, would also be similar to the eigenvalues of R. Particularly, some 
of the extreme eigenvalues of T, would be approximately the same as the 
corresponding eigenvalues of R. Generally, because of the numerical properties of 
the block Lanczos algorithm, it is not practical to carry the method through to 
completion [Ref. 1]. The importance of the algorithm lies in the fact that some of 
the smallest (and largest) eigenvalues of T, will closely approximate the 
corresponding eigenvalues of R for values of s such that qs «« n. This is stated by 
the following theorem. 

Theorem 3.1 [Ref. 1. Let À} < Àp < ++. < À, be the eigenvalues of R and 
Vi, Vo, ***, V, be the corresponding orthonormal eigenvectors. Assume that 
A¿<Aq+1- Apply the block Lanczos recursion in Eqns(3—1)—(3—3) to R generating s 


blocks and let 44 € jt < -+ +< up be the eigenvalues of T,. Suppose that 


W 
w= | VO (3—9) 


Wa 





is the nxq matrix of projections of the starting block of vectors Q, on the 
eigenvectors of R, where V = [v,, Vy, V3, ‘°° , V,| and W, is a qxq matrix composed 


of the first q rows of W. Suppose further that W, is nonsingular so that o the 


min? 
smallest singular value of W,, is greater than zero. Then for k = 1, 2, 3, ..., q, the 
eigenvalues of T; satisfy | 


Ay € y € A. 6 (3-10) 
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where the spread ef is given by 





1 (A, -A,) t an*0 
)وه‎ ) 


and 0 = cos Tonin Yk = (Ak = Aq), ~ Àn) and $,- is the (s—1)"" Chebyshev 
polynomial of the first kind. 

This theorem illustrates the importance of the local gaps |A, — Ags;|, but 
does not show the potential positive effect of the outer loop iteration of an iterative 
block Lanczos procedure on reducing the overall effective spread and thereby 
improving the convergence rate. 

Note that, as we have defined it, the block Lanczos method is not a method for 
finding the eigenvalues and the eigenvectors of a symmetric matrix R. Rather, it is 
a procedure for finding a block tridiagonal matrix T, which is similar to R. To 
produce a complete algorithm for finding the eigenvalues and the eigenvectors, we 
need to combine the Lanczos algorithm with a technique for finding the eigenvalues 
۸ and the eigenvectors y, of T, such as the QR algorithm. 

Now, we will consider certain properties of the block Lanczos algorithm and 
problems associated with its implementation and application. 

The computed Lanczos blocks Q; are desired to be mutually orthogonal. In 
practice, however, because of the arithmetic errors when P; is computed, they 
rapidly lose orthogonality. Thus, after a few iterations of the block Lanczos 
algorithm, the current Q; is no longer orthogonal to the previous Lanczos blocks 


Q,, Qo, ‘°°, Q;- The subsequent losses in orthogonality between the blocks caused 
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by the roundoff errors increase aS we increase the number of Lanczos blocks. At 
some stage when these losses have accumulated sufficiently, the assumption that 
the block tridiagonal Lanczos matrix is the projection of the given matrix R on the 
subspaces Q will be false [Ref. 1]. As a result the q smallest eigenvalues of T, may 
not approximate the q smallest eigenvalues of R. Thus, it requires costly 
reorthogonalization of each Q;,, with respect to all the previous Lanczos blocks to 
maintain the stability of the algorithm [Ref. 4]. 

Loss of orthogonality goes hand-in-hand with the convergence of some of the 
eigenvalues of T, to the eigenvalues of R. In this case we have two options, stop 
the Lanczos iterations as the blocks begin to lose their mutual orthogonality or 
reorthogonalize the blocks if more iterations are desired. The difficulty in using the 
Lanczos method in this way lies in reliably and efficiently determining at what point 
the orthogonality is being lost. In order to reorthogonalize all Lanczos blocks, we 
first reorthogonalize the residual matrix P; with respect to all the previous Lanczos 
blocks and compute the next block Q.,, and the coefficient matrix B,;,, such that 
تل اا‎ This modification preserves the stability of the algorithm but at a 
considerable cost because the reorthogonalization process requires a large number of 
arithmetic operations. The need for reorthogonalization seems to increase with the 
size n of R and the number of blocks that are required to be computed in the 
algorithm [Ref. 1]. 

We now describe how to choose the block size q. It is usually best to choose q 
equal to the number of eigenvalues and eigenvectors r that we are attempting to 
compute which could be the multiplicity of the smallest eigenvalue. Theorem 3.1 
suggests that a good choice for q is one for which the gap between A, and A,,, is 


fairly large. 
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B. ALGORITHM 

It is possible to find a few of the extreme eigenvalues and corresponding 
eigenvectors of a real symmetric matrix using the block Lanczos algorithm rather 
than computing the entire matrix decomposition. Each of these smallest eigenvalues 
and the corresponding eigenvector of the autocorrelation matrix for received signals 
from a sensor array has the spectral information to estimate the 
direction-of-arrival. 

The computer simulation experimental set up used is the same as the one shown 
in Figure 1. We receive the signals at a linear array containing M equally-spaced 


sensors and generate the nxn autocorrelation matrix of these received signals. 


1. Reduction 

The reduction of the data proceeds as follows. Using the block Lanczos 
algorithm we reduce the autocorrelation matrix R into a block tridiagonal matrix 
that has the same extreme eigenvalues as R. In this section we will present a 
practical algorithm to implement the block Lanczos method. 

Given an nxn autocorrelation matrix R, we generate an initial nxq matrix 
Q, which is chosen arbitrarily and orthonormalized. The number of vectors in each 
block, q, is considered to be between 3 and 5 in this study. To begin, we compute 


RQ, and a residual matrix P, given by 


P, = RQ, لا‎ (3-15) 


where A, is a qxq coefficient matrix chosen so that the Euclidean norm of P, is 


minimized with respect to all possible choices of A, [Ref. 4]. It can be shown [Ref. 1] 
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that ||P,l| is minimized when 
1 
Ay = Q;RQ,. (3—16) 


The qxq matrix A, forms the first block on the main-diagonal of the block 


tridiagonal matrix, T, (see Eqn(3—6)). With this choice for A}, we have 


P, = )1 Q,Q})RQ,. (3-17) 


That is, P, is the projection of RQ, onto the space orthogonal to that spanned by 
the columns of Q,. The second Lanczos block of vectors Q, and a qxq upper 
triangular coefficient matrix B, are then obtained by using the QR factorization 


with modified Gram-Schmidt procedure on P;: 


The current Lanczos block Q, is orthonormal to the previous block Q,. The upper 
triangular matrix B, and its transposed version Ba form the first elements in the 
sub-diagonal and the super-diagonal, respectively, of the block tridiagonal matrix 
T 

The remaining matrices in the sequence of the Lanczos blocks وا پا‎ 


Q., where s<<n, are computed as follows: For i=2,3,---,s, compute 


A, 7 QI(RQ; - Q,;BT) (3-19) 


31 


P; = RQ, — Q;A; — Q;.,B} (3-20) 
Q341B;4, = P; (3—21) 


where Q;,, and B;,, are obtained as the QR factorization of the residual matrix P;. 
A modified Gram-Schmidt procedure can be used to reorthgonalize the columns of 
P; This means that Q;,, is orthonormal to all previous matrices Q,, Q5, - - -, Q;. 
As we increase the number of the Lanczos blocks, their mutual orthogonality is 
preserved because of the built-in QR algorithm to factor the residual matrix P. 
Consequently, the space spanned by Q,, Q,, - - -, Q, contains the columns of the 
matrices Q,, RQ,, R?Q,, + + -, R&!'Q, which form the orthonormal basis of the 
Krylov subspace. Thus, we do not need to reorthogonalize the Lanczos blocks when 


we use the algorithm in Eqns(3—17)-(3—19). 


2. | Figendecomposition 

We need to compute the eigenvalues of the block tridiagonal matrix T, 
which approximate the eigenvalues of an autocorrelation matrix. Then the 
eigenvectors of R corresponding to these eigenvalues are determined by knowing the 
matrix of Lanczos blocks. 

There are several techniques for computing the eigenvalues and 
eigenvectors of a given matrix. The fundamental algebraic eigenproblem is to 
determine the eigenvalues 4; given the set of qs homogeneous linear equations in qs 


unknowns [Ref. 17] 


Tey, = Wy;  fori=1,2,---,qS (3-22) 
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where the eigenvectors y; of T, satisfy 


0, Mi], 
1 
A A (3-23) 
. A: 
that is, y; are mutually orthonormal vectors. From Eqn(3—22) the characteristic 


equation associated with the matrix T, is given by 
det(T, — 1,1) = 0. (3-24) 
Expanding the determinant, we have the polynomial equation 
A O AO (3-25) 


where a; are the coefficients of jj and the roots of this polynomial give us the 
eigenvalues jj of the T.. Corresponding to any eigenvalue 4;  Eqn(3—22) has at 
least one non-trivial solution y;. Since the eigenvalues of T, approximate a few of 
the extreme eigenvalues of autocorrelation matrix R, we choose the q smallest 


eigenvalues and the corresponding eigenvectors of T, and use the Ritz vector to 


compute the eigenvectors of R given by [Ref. 1] 

X = QY (3-26) 
where X is a group of approximate eigenvectors of R, X=[x, x) --- xj, Q is a 
matrix of Lanczos blocks, and Y is a group of q smallest eigenvectors of T,, Y=[y, 


Yo °° ا‎ Also, we can use the Rayleigh quotient iteration to find the eigenvectors 
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of R related to the eigenvalues of T,. However, this method is slower than using 
the Ritz vector because it needs more time to converge for a given eigenvalue. 

Now, using the eigenvectors of the autocorrelation matrix R we estimate 
the spatial spectrum of the received signals. Each of the eigenvectors contains the 
true spectral information as well as some spurious peaks. The direction-of-arrival 
of point sources can be estimated by computing the spatial power spectral density of 
the eigenvectors of R. The power spectral density estimate for the 7" eigenvector 


corresponding to one of the smallest eigenvalues of R is computed as follows: 


SON M (3-27) 
وګ‎ | es 


where x;; are the elements of the j'® eigenvector, x; and 0 < f< 0.5. 
C. SIMULATION RESULTS 

Using the block Lanczos algorithm we can selectively compute a few of the 
smallest eigenvalues and eigenvectors of an autocorrelation matrix. These 
eigenpairs are in the noise subspace and contain the spectral information of the 
source bearings from an array of sensors. Thus, we could estimate the spatial 
power spectral density for each eigenvector using Eqn(3—27). We have used the 
spectral product of several individual PSDs to improve the direction-of-arrival 
estimation performance. The advantage in using the multiplicative PSD function 
was demonstrated in Chapter 2 for the single vector Lanczos algorithm. The same 


advantage holds for the block Lanczos algorithm. 
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The computer simulation experiment consists of an equally-spaced linear array 
of 100 sensors receiving signals of known temporal frequency from various bearings. 
The size of the autocorrelation matrix computed is 25x25. The number of smallest 
eigenpairs to be estimated is q, which is the size of the Lanczos block. We choose q 
to be 3, 5, or 7 in these examples and used 5 dB, 0 dB, —5 dB and —10 dB 
signal-to-noise ratios (SNR). 

Example 1 is the detection of a target at 9° for different SNRs. Figure 10 shows 
the spectral overlay of 5 eigenvectors corresponding to the smallest eigenvalues at 
an SNR of 5 dB. The spectrum of each eigenvector has a common peak at the true 
bearing 9°. Figure 11(a) and Figure 11(b) show improved DOA estimation where 
the former figure illustrates the average of 5 eigenvectors and the latter figure shows 
the spectral multiplication of those eigenvectors. As can be seen, the spectral 
multiplication technique has greatly improved the nulls and the peak at 9°. Thus, 
in the remainder of the results, we have used the spectral multiplication method to 
compute the PSD function. Figure 12(a) shows the performance at 0 dB using 3 
eigenvectors and Figure 12(b) shows the result at an SNR of 0 dB using 5 
eigenvectors. Improved performance is obtained by using more eigenvectors. In 
Figure 12(a) we have several large spurious peaks around the true peak at 9°, but 
in Figure 12(b) we have just 2 small spurious peaks and the true bearing is clearly 
evident. Figure 13(a) and Figure 13(b) show the results at an SNR of —5 dB using 3 
and 5 eigenvectors, respectively. The spurious peaks have increased magnitudes 
but the difference between the true peak and the spurious peaks is still large enough 
to determine the true bearing. Figure 14(a) and Figure 14(b) show the result at an 
SNR of —10 dB. Using 3 eigenvectors (Figure 14(a)), it is hard to determine the 


true bearing. Using 5 eigenvectors (Figure 14(b)), the true peak appears to be 


4 


larger than the spurious peaks and we can recognize the true bearing. 

In Example 2 we have 3 targets at 34°, 36° and 54°. Two targets are very 
closely spaced in bearing. At 0 dB, the results of two cases, where one has used 5 
eigenvectors (Figure 15(a)) and the other 7 eigenvectors (Figure 15(b)), indicate 
good performance. There are no spurious peaks and the resolution is clearly 
acheived. Figure 16(a) shows the result at —5 dB using 5 eigenvectors. In this 
result, the resolution is acheived but one spurious peak is as large as the true peaks. 
In Figure 16(b), the result shows improved performance using 7 eigenvectors. 
There are no spurious peaks and good spectral resolution is acheived. At —10 dB, 
the algorithm cannot separate the two targets located at 34° and 36° when we use 5 
eigenvectors (Figure 17(a)). When we use 7 eigenvectors however (Figure 17(b)) it 
can separate the two closely located sources. Nevertheless, a number of spurious 
peaks are larger in magnitude than the peak at 36° making it impossible to 
determine the true bearings accurately. 

The results in this chapter show that the eigenvectors computed using the block 
Lanczos algorithm can be used to determine the spectrum even in very low SNRs. 
Since the block Lanczos method can compute a few of the extreme eigenpairs of a 


large symmetric matrix, it is efficient in computations and storage. 
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D. COMPARISON 

In this section we compare the performance of the block Lanczos algorithm and 
the single vector Lanczos algorithm for DOA estimation. In both cases, we chose 
the seven smallest eigenvalues and corresponding eigenvectors and use the spectral 
multiplication method since this gave the best preformance for both algorithms. 

Example 1 is the comparison for detection of two targets at 18° and 45° for 
different SNRs. Figure 18(a) shows the result using the single vector Lanczos 
algorithm and Figure 18(b) indicates the result using the block Lanczos algorithm at 
an SNR of 0 dB. The result of the single vector case has two small spurious peaks; 
the block case has one very small spurious peak. In both cases the true bearings are 
clearly distinguished from the spurious peaks. Figure 19 shows a comparison of the 
results at —5 dB. The performance of the block Lanczos algorithm (Figure 19(b)) is 
much better than the single vector case (Figure 19(a)) even though several large 
spurious peaks are present. Figure 20(a) is the result at —10 dB using the single 
vector Lanczos algorithm and Figure 20(b) is the corresponding result of the block 
case. In the single vector case, the spurious peaks are almost the same as the true 
peaks and we cannot distinguish the true bearing location. In the block case, 
however, the true peaks are slightly larger than the largest spurious peak. 

Example 2 is a comparison for detection of four targets at 18°, 27°, 29°, and 
45°. Figure 21 shows the performance of the two methods at 0 dB. The resolution 
is clearly achieved in both cases. Figure 22 shows the results at —5 dB. 
Figure 22(a) shows that the single vector Lanczos algorithm cannot separate the two 
closly spaced targets. The block case (Figure 22(b)) shows slightly better resolution 
between targets located at 27° and 29°. At —10 dB (Figure 23(a) and Figure 23(b)) 


neither method can separate the two closely spaced in bearings at 27° and 29°. 
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Also, a spurious peak is produced by both methods. 

The block Lanczos algorithm is known to estimate the multiple eigenvalues 
better than the single vector case [Ref. 1]. This situation is applicable to the DOA 
estimation using noise subspace computation where the noise is white. As observed 
in the results of Figures 18 — 23, the spectral estimation performance of the block 
Lanczos algorithm is consistently better than the single vector algorithm, 
particulary at low SNRs. Also, the block method provided better spectral 
resolution than the single vector method in our tests. While more analysis is needed 
to validate the simulation results, the overall performance of the block method is 


quite encouraging. 
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Figure 18. "Two targets at 18° and 45° (0 dB) : 


(a) single vector Lanczos method (b) block Lanczos method 
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Figure 19. Two targets at 18° and 45° (-5 dB) : 


(a) single vector Lanczos method (b) block Lanczos method 


94 


10° = 


10° 


107 


10 


105 


Power Spectral Density 


Bearing 


Power Spectal Density 





Bearings 


Figure 20. Two targets at 18° and 45° (-10 dB) : 


(a) single vector Lanczos method (b) block Lanczos method 
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Figure 21. Four targets at 18°, 27°, 29° and 45° (0 dB) : 


(a) single vector Lanczos method (b) block Lanczos method 
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Figure 22. Four targets at 18°, 27°, 29° and 45° (-5 dB) : 


(a) single vector Lanczos method (b) block Lanczos method 
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Figure 23. Four targets at 18%, 27%, 29% and 45% (-10 dB) : 


(a) single vector Lanczos method (b) block Lanczos method 
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IV. SUMMARY AND CONCLUSIONS 


In this thesis we examined the use of the single vector Lanczos method and the 
block Lanczos method and its application to spectral analysis and 
direction-of-arrival problems. 

We computed a few of the extreme eigenvalues and their associated 
eigenvectors of a large symmetric matrix using the block Lanczos method. The 


eigenvalues and eigenvectors of the Lanczos matrix T., approximate the 


5 
corresponding eigenvalues and eigenvectors of the given matrix R. The block 
Lanczos algorithm can directly determine the multiplicities of the effective 
eigenvalues and the eigenvectors of R. We found that the spectral estimate of the 
block Lanczos method is more accurate than the single vector Lanczos method, 
particulary at low SNRs. Since we compute only a few of the extreme eigenpairs of 
a large autocorrelation matrix, the result of this algorithm is savings in 
computations and storage. This algorithm may be applied to any system where one 
needs to obtain a rapid decomposition of a large correlation matrix. 

Although the results of this thesis are most encouraging, some additional 
work still remains to be done. We need to compare the results in computational 
speed and accuracy with other eigendecomposition techniques for valdating this 


algorithm. Also, we need to analyze the effect of roundoff errors on the eigenpair 


estimation of the Lanczos algorithm. 
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